ICPC Live Archive/3928 - Ballroom lights/help/gomox.ar-pap-2010-d9b17e5cb110/tpd/ballroom/ballroom.py

   1 #!/usr/bin/python
   2 # -*- coding: utf8 -*-
   3
   4 from __future__ import division
   5
   6 import psyco
   7 psyco.full()
   8
   9 import math
  10
  11
  12 ###########
  13 # Helpers #
  14 ###########
  15
  16 def mod2pi(x):
  17     '''Devuelve un valor de angulo equivalente en [0, 2*pi).
  18
  19     >>> import math
  20     >>> mod2pi(0.0)
  21     0.0
  22     >>> mod2pi(4 * math.pi)
  23     0.0
  24     >>> mod2pi(-math.pi)
  25     3.1415926535897931
  26     '''
  27     while x < 0:
  28          x += 2 * math.pi
  29     while x >= 2 * math.pi:
  30         x -= 2 * math.pi
  31     return x
  32
  33
  34 def mod360(x):
  35     '''Devuelve un valor de angulo equivalente en [0, 360).
  36
  37     >>> import math
  38     >>> mod360(0.0)
  39     0.0
  40     >>> mod360(360+180)
  41     180.0
  42     >>> mod360(-180)
  43     180.0
  44     '''
  45     while x < 0:
  46          x += 360.0
  47     while x >= 360.0 :
  48         x -= 360.0
  49     return x
  50
  51
  52 ##########
  53 # Vector #
  54 ##########
  55
  56 class Vector:
  57     def __init__(self, x, y):
  58         self.x = x
  59         self.y = y
  60
  61     def __add__(self, otro):
  62         return Vector(self.x + otro.x, self.y + otro.y)
  63
  64     def __sub__(self, otro):
  65         return Vector(self.x - otro.x, self.y - otro.y)
  66
  67     def __mul__(self, f):
  68         return Vector(self.x * f, self.y * f)
  69
  70     def __rmul__(self, f):
  71         return self * f
  72
  73     def __div__(self, f):
  74         return Vector(self.x / f, self.y / f)
  75
  76     def __neg__(self):
  77         return Vector(self.x *-1 , self.y * -1)
  78
  79     def __eq__(self, otro):
  80         return self.x == otro.x and self.y == otro.y
  81
  82     def __repr__(self):
  83         return '''<Vector %s,%s>''' % (self.x, self.y)
  84
  85     def __iter__(self):
  86         yield self.x
  87         yield self.y
  88
  89     def __hash__(self):
  90         return hash((self.x,self.y))
  91
  92     def dotP(self, otro):
  93         '''Computa el producto escalar de 2 vectores.
  94
  95         >>> Vector(1,0).dotP(Vector(0,1))
  96         0.0
  97         >>> Vector(1,1).dotP(Vector(1,1))
  98         2.0
  99         '''
 100
 101         return self.x * otro.x + self.y * otro.y
 102
 103     def boundNorm(self, norma):
 104         '''Devuelve un vector idéntico a este pero con la norma
 105         acotada al rango [0, norma]
 106
 107         >>> Vector(2,0).boundNorm(1)
 108         <Vector 1.0,0.0>
 109         >>> Vector(1,1).boundNorm(4)
 110         <Vector 1.0,1.0>
 111         >>> Vector(2,2).boundNorm(math.sqrt(2))
 112         <Vector 1.0,1.0>
 113         '''
 114
 115         assert norma > 0
 116         norma_efectiva = min(norma, self.norma())
 117         factor = norma_efectiva / self.norma()
 118         return self * factor
 119
 120     def flipY(self):
 121         '''Devuelve un vector idéntico a este con la coordenada
 122         Y invertida (útil para dibujar en la pantalla, que tiene el
 123         sistema de coordenadas con esta particularidad).
 124         '''
 125         return Vector(self.x, -self.y)
 126
 127     def norma(self):
 128         '''Computa la norma cartesiana del vector.
 129
 130         >>> Vector(1,1).norma()
 131         1.4142135623730951
 132         >>> Vector(1,0).norma()
 133         1.0
 134         >>> Vector(0,1).norma()
 135         1.0
 136         '''
 137         return math.sqrt(self.x**2 + self.y**2)
 138
 139     def angulo(self):
 140         '''Computa el ángulo polar del vector (en radianes desde la recta y=0).
 141
 142         >>> Vector(1,0).angulo()
 143         0.0
 144         >>> Vector(0,1).angulo()
 145         1.5707963267948966
 146         >>> Vector(0,-1).angulo()
 147         -1.5707963267948966
 148         '''
 149         if self.x == 0 and self.y == 0:
 150             return 0
 151         elif self.x >= 0:
 152             return math.asin(self.y/self.norma())
 153         else:
 154             return -math.asin(self.y/self.norma()) + math.pi
 155
 156     def grados(self):
 157         '''Computa el ángulo polar del vector (en grados desde la recta y=0).
 158
 159         >>> Vector(1,0).grados()
 160         0.0
 161         >>> Vector(0,1).grados()
 162         90.0
 163         >>> Vector(0,-1).grados()
 164         -90.0
 165         '''
 166         return math.degrees(self.angulo())
 167
 168     def viento(self):
 169         '''Devuelve el viento correspondiente a este vector.
 170
 171         >>> viento = Vector(10,10).viento()
 172         >>> viento
 173         <Viento 225.0º de 14.1kt>
 174         >>> Vector(1, 0).viento()
 175         <Viento 270.0º de 1.0kt>
 176         >>> Vector(-10, 0).viento()
 177         <Viento 90.0º de 10.0kt>
 178         '''
 179         return Viento(mod360(270 - self.grados()), self.norma())
 180
 181
 182 ##########
 183 # Viento #
 184 ##########
 185
 186 class Viento:
 187     def __init__(self, direccion, intensidad):
 188         '''Construye un viento soplando desde la direccion indicada,
 189         con la intensidad indicada.
 190
 191         >>> viento = Viento(0, 15)
 192         >>> viento.angulo
 193         4.7123889803846897
 194         >>> viento.x
 195         -2.7553642961003488e-15
 196         >>> viento.y
 197         -15.0
 198         '''
 199
 200         assert 0 <= intensidad
 201         direccion = mod360(direccion)
 202         self.direccion = direccion
 203         self.intensidad = intensidad
 204
 205         # Precalculo algunas cosas
 206         self.angulo = mod2pi(-math.radians(direccion) - math.pi/2)
 207         self.x = math.cos(self.angulo) * intensidad
 208         self.y = math.sin(self.angulo) * intensidad
 209
 210     def vector(self):
 211         '''Devuelve el vector asociado a este viento.
 212
 213         >>> vector = Viento(90,10).vector()
 214         >>> vector.x
 215         -10.0
 216         >>> vector.y
 217         1.2246063538223773e-15
 218         >>> vector.viento()
 219         <Viento 90.0º de 10.0kt>
 220         '''
 221         return Vector(self.x, self.y)
 222
 223     def anguloIncidencia(self, vector_navegacion):
 224         '''Devuelve el ángulo de incidencia del viento sobre un
 225         barco que navega sobre vector_navegacion.
 226
 227         >>> Viento(0, 10).anguloIncidencia(Vector(-1,1))
 228         0.78539816339744828
 229         >>> Viento(0, 10).anguloIncidencia(Vector(1,1))
 230         5.497787143782138
 231         '''
 232         return mod2pi(vector_navegacion.angulo() - (-self.vector()).angulo())
 233
 234     def __repr__(self):
 235         return "<Viento %.1fº de %.1fkt>" % (self.direccion, self.intensidad)
 236
 237
 238
 239
 240
 241 inf = float('inf')
 242
 243 def intersect_lines(l1, l2):
 244
 245     p1 = l1[0]
 246     p2 = l1[1]
 247     p3 = l2[0]
 248     p4 = l2[1]
 249
 250     den = (p4.y-p3.y)*(p2.x-p1.x) - (p4.x-p3.x)*(p2.y-p1.y)
 251
 252     ua_num = (p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x)
 253     ub_num = (p2.x-p1.x)*(p1.y-p3.y) - (p2.y-p1.y)*(p1.x-p3.x)
 254
 255     if abs(den) == 0:
 256         if abs(ua_num - ub_num) == 0 and \
 257            abs(ub_num) == 0:
 258             return max(p1, p2, key=lambda v: v.norma())
 259
 260         return None
 261
 262     ua = ua_num/den
 263     ub = ub_num/den
 264
 265     if ub >= 0:
 266         return ua
 267
 268     return None
 269
 270
 271
 272 def rango_oscurecido(inicio, fin, src, centro, radio):
 273     '''
 274     Devuelve el rango (en distancia de 0 a ||$inicio-$fin||)
 275     que se encuentra a la sombra de una columna de radio $radio,
 276     ubicada en $centro e iluminada por una fuente de luz puntual en $src.
 277     '''
 278
 279     bisectriz = centro - src
 280     p = math.asin(radio/bisectriz.norma())
 281     alpha = bisectriz.angulo()
 282
 283     d1 = Vector(math.cos(alpha + p), math.sin(alpha + p))
 284     d2 = Vector(math.cos(alpha - p), math.sin(alpha - p))
 285
 286     p1 = intersect_lines((inicio, fin), (src, src + d1))
 287     p2 = intersect_lines((inicio, fin), (src, src + d2))
 288
 289     unidad = (fin - inicio).norma()
 290
 291     # Se sabe que los puntos de interseccion estan ordenados correctamente
 292     # porque el angulo de diferencia entre las semirectas es menor a 180,
 293     # entonces no es necesario chequearlo.
 294     if p1 is None:
 295         if p2 is None:
 296             return (0.0, 0.0)
 297         else:
 298             res = (0.0, p2)
 299     else:
 300         if p2 is None:
 301             res = (p1, 1.0)
 302         else:
 303             res = (p1, p2)
 304
 305     res_bounded = interseccion_intervalos(res, (0.0, 1.0))
 306     p = (res_bounded[0] * unidad, res_bounded[1] * unidad)
 307     return p
 308
 309
 310 def interseccion_intervalos(i1, i2):
 311     '''Devuelve los extremos del intervalo interseccion entre i1 y i2.
 312
 313     >>> interseccion_intervalos((1,3),(2,8))
 314     (2, 3)
 315     >>> interseccion_intervalos((1,2), (5,6))
 316     (5, 5)
 317     >>> interseccion_intervalos((-inf, 8), (6, inf))
 318     (6, 8)
 319     '''
 320
 321     p1, p2 = i1
 322     p3, p4 = i2
 323
 324     d = (max(p1, p3), min(p2, p4))
 325     if d[1] < d[0]:
 326         return (d[0], d[0])
 327     else:
 328         return d
 329
 330
 331 def reduce_interseccion(l):
 332     cand = (-inf, inf)
 333     for e in l:
 334         cand = interseccion_intervalos(cand, e)
 335     return cand
 336
 337
 338 def reduce_union(l):
 339     '''Devuelve la union de la lista de intervalos l
 340
 341     >>> reduce_union([(1,3), (2,5)])
 342     [(1, 5)]
 343     >>> reduce_union([(1,3), (4, 6)])
 344     [(1, 3), (4, 6)]
 345     >>> reduce_union([(1,3), (4, 6), (5, 6)])
 346     [(1, 3), (4, 6)]
 347     >>> reduce_union([(1,3), (4, 6), (5, 7)])
 348     [(1, 3), (4, 7)]
 349     '''
 350
 351     ordenados = l[:]
 352     ordenados.sort(key=lambda x: x[0])
 353
 354     if len(l) == 0:
 355         return []
 356
 357     res = []
 358     cand = ordenados[0]
 359     for i in range(1, len(l)):
 360         otro = ordenados[i]
 361         if otro[0] > cand[1] and otro[0] != otro[1]:
 362             res.append(cand)
 363             cand = otro
 364         else:
 365             if otro[1] > cand[1]:
 366                 cand = (cand[0], otro[1])
 367
 368     res.append(cand)
 369     return res
 370
 371
 372 # Asume que los intervalos estan ordenados por punto de comienzo
 373 def invertir_rango(ints, fin):
 374     '''
 375     Invierte el conjunto de segmentos (minimo y ordenado) para
 376     tomar su complemento.
 377
 378     >>> invertir_rango([(3,4)], 8)
 379     [(0, 3), (4, 8)]
 380     >>> invertir_rango([(3,4),(6,8)], 8)
 381     [(0, 3), (4, 6)]
 382     >>> invertir_rango([(0,2), (5,6)], 10)
 383     [(2, 5), (6, 10)]
 384     '''
 385
 386     if len(ints) == 0:
 387         return [(0, fin)]
 388
 389     invertido = []
 390     if ints[0][0] > 0:
 391         invertido.append((0, ints[0][0]))
 392
 393     for i in range(1, len(ints)):
 394         invertido.append((ints[i-1][1], ints[i][0]))
 395
 396     if ints[-1][1] < fin:
 397         invertido.append((ints[-1][1], fin))
 398
 399     return invertido
 400
 401
 402
 403
 404
 405 def resolver(lamparas, columnas, max_x, max_y):
 406
 407     paredes = ((Vector(0,0), Vector(0,max_y)),
 408                (Vector(0, max_y), Vector(max_x, max_y)),
 409                (Vector(max_x, max_y), Vector(max_x, 0)),
 410                (Vector(max_x, 0), Vector(0,0)))
 411
 412
 413     total = 0
 414     for p in paredes:
 415         luces_pared = []
 416         for l in lamparas:
 417             zonas_tapadas = []
 418             for c in columnas:
 419                 r = rango_oscurecido(p[0], p[1], l, c[0], c[1])
 420                 zonas_tapadas.append(r)
 421             zona_iluminada = invertir_rango(reduce_union(zonas_tapadas), (p[1]-p[0]).norma())
 422             luces_pared += zona_iluminada
 423         luces_pared = reduce_union(luces_pared)
 424
 425         s = sum(map(lambda x: x[1] - x[0], luces_pared))
 426         total += s
 427
 428     print "%.4f" % total
 429
 430
 431
 432 def main():
 433     import sys
 434
 435     while True:
 436         nl, nc, max_x, max_y = map(int, sys.stdin.readline().strip().split())
 437
 438         if nl == nc == max_x == max_y == 0:
 439             break
 440
 441         lamparas = []
 442         for i in range(nl):
 443             lampara = Vector(*(map(float, sys.stdin.readline().strip().split())))
 444             lamparas.append(lampara)
 445
 446         columnas = []
 447         for i in range(nc):
 448             x, y, r = map(float, sys.stdin.readline().strip().split())
 449             columnas.append((Vector(x,y), r))
 450
 451         resolver(lamparas, columnas, max_x, max_y)
 452
 453
 454
 455 if __name__ == '__main__':
 456     #import doctest
 457     #doctest.testmod()
 458
 459     main()